% Calculate temperature sensitivity from shot noise expectation and compare
% to experiment
% hyzhou, 11/22/2017
% Input:
%   FullListfl: full fluorescence counter data
%   gSeqParameter: certain acquistion parameters for the fluorescence
%   MeasTime: measurement times
% Output:
%   theor: shot noise limit
%   exp: experimental sensitivity
%   relfluct: relative fluctuations normalized by shot noise expectation
function [theor, exp, relfluct] = get_shot_noise_limit(FullListfl,gSeqParameter,MeasTime, varargin)

bTogether = 0;
bPlot = 0;
fontsize = 16;
if ~isempty(varargin)
    bTogether = varargin{1}(1);
    bPlot = varargin{1}(2);
    fontsize = varargin{1}(3);
end
% calculate acquisition time for each of these data counters
% check whether differences between the counters within 1 cycle is
% consistent with shot noise
int_time = gSeqParameter.NSamples * gSeqParameter.GreenOn * 1e-9;   % ns
fl = FullListfl*int_time*1e6;   % Mcps
if size(fl,1) > 4
    fldiffs = reshape(fl(1:4,:) - fl(8:-1:5,:), [4*size(fl,2),1]);
    flaves = reshape(fl(1:4,:) + fl(8:-1:5,:), [4*size(fl,2),1]);
    relfluct = fldiffs./sqrt(flaves);
    fl = (FullListfl(1:4,:) + FullListfl(8:-1:5,:))/2;
else
    relfluct = 0;
end

% Theoretically expected sensitivity
cyc_time = MeasTime(end)/length(MeasTime);
% green is on twice
count_rate = gSeqParameter.NSamples * gSeqParameter.GreenOn * 2 * 1e-9 * mean(mean(FullListfl)) * gSeqParameter.reps * (gSeqParameter.Navrg/2) * 1e6;
theor = sqrt(2) * sqrt(cyc_time) * gSeqParameter.delta_w / (gSeqParameter.sens_fin*2 * sqrt(count_rate) * 74e3);

% Experiment duty cycle for acquisition, 2 is for IR on/off, 4 is for MW
% counters
%duty_cyc = count_rate / mean(mean(FullListfl*1e6)) * 2 * 4 / cyc_time;

% Experimentally extracted sensitivity
if bTogether
    % put same reps together
    fl1 = reshape(fl,[4,gSeqParameter.reps,size(fl,2)/gSeqParameter.reps]);
    fl2 = reshape(mean(fl1, 2),[4,size(fl,2)/gSeqParameter.reps]);
    fl = fl2;
else
    % alternatively, separate each repetition
    reps = gSeqParameter.reps;
    fl1 = reshape(fl,[4,4*reps,size(fl,2)/reps/4]);
    temp = [(1:reps)',(2*reps:-1:reps+1)',(2*reps+1:3*reps)',(4*reps:-1:3*reps+1)'];
    temp = reshape(temp',[4*reps,1]);
    fl2 = fl1(:,temp,:);
    fl = reshape(fl2,[4,size(fl,2)]);
end

ind_list = 1:size(fl,2);
ind_on = ind_list(mod(ind_list,4)==2 | mod(ind_list,4)==3);
ind_off = ind_list(mod(ind_list,4)==0 | mod(ind_list,4)==1);
fl_on = fl(:,ind_on);
fl_off = fl(:,ind_off);
bin_vals = round(logspace(0,1.4,10))'; % adjust here depending on number of temperature data: round(logspace(0,2,10))
std_vals = zeros(size(bin_vals));
temp_vals = zeros(size(bin_vals));
for j = 1:length(bin_vals)
    bin_size = bin_vals(j);
    fl_on_t = fl_on(:,1:(size(fl_on,2)-mod(size(fl_on,2),bin_size)));
    fl_on_t = mean(reshape(fl_on_t, [4,size(fl_on_t,2)/bin_size,bin_size]),3);
    temp_on = ((fl_on_t(1,:) + fl_on_t(2,:)) - (fl_on_t(3,:)+ fl_on_t(4,:))) ./ ((fl_on_t(1,:) - fl_on_t(2,:)) + (fl_on_t(4,:)- fl_on_t(3,:)))*gSeqParameter.delta_w/(-gSeqParameter.t_slope);
    fl_off_t = fl_off(:,1:(size(fl_off,2)-mod(size(fl_off,2),bin_size)));
    fl_off_t = mean(reshape(fl_off_t, [4,size(fl_off_t,2)/bin_size,bin_size]),3);
    temp_off = ((fl_off_t(1,:) + fl_off_t(2,:)) - (fl_off_t(3,:)+ fl_off_t(4,:))) ./ ((fl_off_t(1,:) - fl_off_t(2,:)) + (fl_off_t(4,:)- fl_off_t(3,:)))*gSeqParameter.delta_w/(-gSeqParameter.t_slope);
    std_vals(j) = std(temp_on-temp_off);
    temp_vals(j) = mean(temp_on-temp_off);
end
ft1 = fittype( 'a*x^(-0.5)', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
[fitresult, ~] = fit(bin_vals*MeasTime(end)/length(fl_on), std_vals, ft1, opts);
exp = fitresult.a;

if bPlot
    figure(201)
    set(gcf,'position',[100,100,800,800])
    clf
    subplot(2,1,1)
    dstep = 0.1;
    xx = -5:dstep:5;
    hist(relfluct,xx)
    hold on
    plot(xx,length(relfluct)*dstep*normpdf(xx,0,1))
    xlabel('Relative difference', 'fontsize', fontsize)
    ylabel('Frequency', 'fontsize', fontsize)
    set(gca,'fontsize', fontsize)
    title(['IR current: ',num2str(gSeqParameter.InfraredmA),' mA, MW power: ',...
        num2str(gSeqParameter.Power),' dBm'])
    subplot(2,1,2)
    plot(bin_vals*MeasTime(end)/length(fl_on),std_vals,'o-')
    hold on
    plot(fitresult, bin_vals*MeasTime(end)/length(fl_on), std_vals);
    title(['Sensitivity: ',...
        num2str(fitresult.a),' K/sqrt(Hz), Theory: ', num2str(theor), ' K/sqrt(Hz)'])
    legend('off')
    xlabel('Integration time', 'fontsize', fontsize)
    ylabel('Standard deviation', 'fontsize', fontsize)
    set(gca,'xscale','log','yscale','log','fontsize', fontsize)
end